Astronomy & Astrophysics manuscript no. shifted convergence zone 


©ESO 2013 


February 13, 2013 





Convergence zones for Type I migration: an inward shift for 

multiple planet systems 

Christophe Cossou^'^, Sean N. Raymond^'^, and Amaud Pierens^'^ 

' Univ. Bordeaux, LAB, UMR 5804, F-33270, Floirac, France. 
~ CNRS, LAB, UMR 5804, F-33270, Floirac, France 

5 December 2012/7 February 2013 

ABSTRACT 

Earth-mass planets embedded in gaseous protoplanetary disks undergo Type I orbital migration. In radiative disks an additional 
component of the corotation torque scaling with the entropy gradient across the horseshoe region can counteract the general inward 
migration. Type I migration can then be directed either inward or outward depending on the local disk properties. Thus, special 
locations exist in the disk toward which planets migrate in a convergent way. Here we present N-body simulations of the convergent 
migration of systems of low-mass (M=l-10Ma) planets. We show that planets do not actually converge in convergence zones. Rather, 
they become trapped in chains of mean motion resonances (MMRs). This causes the planets' eccentricities to increase to high enough 
values to alfect the structure of the horseshoe region and weaken the positive corotation torque. The zero-torque equilibrium point 
of the resonant chain of planets is determined by the sum of the attenuated corotation torques and unattenuated differential Lindblad 
torques acting on each planet. The effective convergence zone is shifted inward. Systems with several planets can experience stochastic 
migration as a whole due to continuous perturbations from planets entering and leaving resonances. 

Key words. Planets and satellites: formation, Protoplanetary disks, Planet-disk interactions, planetary systems. Methods: numerical 



1. Introduction 

An Earth- to Neptune-mass planet embedded in a gaseous pro- 
toplanetary disk excites density waves in the disk (|GoIdreich & 
|Tremaine|1979p , and the back-reaction of these perturbations on 
the planet's orbit causes it to undergo so-called Type I migration 
(Ward 1997[ l. In isothermal disks Type I migration is driven only 
by the differential Lindblad torque and is inward-directed and 
fast (Tanaka et al. 2002| i. In radiative disks, a torque arising from 
material in a planet's horseshoe region can counteract the differ- 
ential Lindblad torque such that migration can be directed either 
inward or outward ( Paardekooper & Mellema; ( ,2006^ ; |Kley & 
|Crida| ( |2008| l). There are locations within disks for which the 
migration is convergent, and we refer to these as convergence 
zones (CZs; Lyra et al.|2010 Paardekooper et al.|201 1\ . 

At a convergence zone, the positive corotation torque is in 
exact balance with the negative differential Lindblad torque, 
such that the total torque experienced by the planet is zero. It 
has been proposed that CZs may act to concentrate planetary em- 
bryos and build large cores ( Lyra et al.|2010 Horn et al.|2012[ l. 
As planets migrate toward the CZ they become trapped in mean- 
motion resonances (MMRs), which prevent unlimited accretion 
( |Morbidelli eraL| (2008 ), S andor et alI] ( [20TT) ). However, colH- 
sions do occur once embryos are trapped in a long enough res- 
onant chain to destabilize it. In addition, turbulence may act to 
break resonances and enhance accretion. 



|Bitsch & Kleyl ( |2010| l show that the corotation torque is at- 
tenuated when a planet has an eccentricity such that its orbit 
crosses a significant fraction of the width of the horseshoe re- 
gion. When two planets become trapped in resonance due to 
convergent migration they excite and sustain non-zero eccen- 
tricities despite eccentricity damping by the disk (e.g. Cresswell 
|& Nelson|2008| l. This in turn should decrease the magnitude of 



the corotation torque and thus modify the balance between the 
differential Lindblad and corotation torques and the subsequent 
migration. 

Here we present simulations of the convergent migration of 
low-mass (M=l-10Me) planets in idealized gas disks. We use a 
simple model for the influence of the eccentricity on the coro- 
tation torque. We show that planets do not actually converge at 
the convergence zone but are instead systematically driven in- 
ward to a shifted equilibrium position with zero net torque. The 
location of the equilibrium position depends on the eccentricity 
excitation of the planets due to mutual perturbations. 

The paper is laid out as follows. In §2 we present our numer- 
ical methods and assumptions. In §3 we present results for two 
planets. In §4 we present results for 3-10 planets. We discuss our 
results in §5. 



2. Methods 

We use an artificial convergence zone (CZ) that mimics a mass- 
independent CZ (whose position does not depend on the mass 
of the planet) at an opacity transition, see (left panel of) Fig.[T] 
where there is an abrupt change from outward to inward migra- 
tion (see, e.g., |Masset|201 1 ). A step function was not used be- 
cause the step in the 'real' profile is only due to the opacity table 
not being smoothed. The location of the CZ was 3 AU. Inside 3 

AU the total torque is positive and equal to Fq = S^rp^Q^,^. 
Outside 3 AU the total torque is -Fq. Here q is the ratio between 
the mass of the planet and the star, h is the aspect ratio that de- 
pend on the temperature profile but is typically 0.05, S^, r^, and 
Qp are respectively the surface density, the orbital distance, and 
the angular speed for the planet. The total torque is a sum of the 
differential Lindblad torque F^ — assumed to remain constant 
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— and the corotation torque Tc- The main interest of the artifi- 
cial CZ is to get rid of the very complex shape of the real profile, 
to only retain the CZ, and nothing more. 

|Bitsch & K ley (2010) show that as the eccentricity of a body 
increases, the structure of its horseshoe region is modified and 
Tc decreases. We implemented a simple formula that reproduces 
the general influence of the eccentricity on Fc by a simple fit to 
the 3D simulations of ( ^Bitsch & Kley|2010| l: 



D = 



rc(g) 

Tc(e = 0) 



- 1 + a- 



tanh(c) - tanh 



lb * e 



+ c 



(1) 



where represents the half-width of the horseshoe region di- 
vided by the semi major axis, e is the planet's eccentricity, and 
our fit produced a - 0.45, b - 3.46, and c - -2.34. We define 
Xs as (see Paardekooper et al.|2010 eq. 44): 



1.1 



r 



1/4 



(2) 



where y is the adiabatic index, q the planet-to-stellar mass ratio, 
and h the disk aspect ratio. The right hand panel of Fig. [T] shows 
that our simple equation is a good match to hydrodynamical sim- 
ulations of the eff'ect of e on Fc, especially at small e, although 
the simulation data points are sparse and appear to include some 
random fluctuations. 
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Fig. 1. Left: The total torque profile of our standard disk in units 
of Fo is shown in red. The dashed blue curve shows the true 
torque profile for a \Qm^ planet near an opacity transition, cal- 
culated using the equations in |Paardekooper et"ar ( 201 1 ). Right: 
Decrease in the corotation torque F^ as a function of planetary 
eccentricity e. We assume that the damping (0 < D < 1) of 
the corotation torque as a function of eccentricity is the same 
for isothermal and fully radiative disc. Then, we derived D from 
Fig. 2 of .Bitsch & Kley ( 2010) by taking the difference between 
their fully radiative and isothermal case, and normalizing to 1 
the "e - case". 



To perform our simulations we used a modified version of 
the Mercury integrator (Chamb ers|1999| ). We included both the 
artificial convergence zone (where F = F^ -i- F^) and F^ varies 
with the planet's eccentricity following Eq{T] Eccentricity and 
inclination damping is done using formulas from Cresswell &| 
jNelso n (2008 ). We assumed the presence of a gas disk with sur- 
face density profile E = 500 (r/1 AU)"'^^ g.cm"^ (used to calcu- 
late Fo and the damping). 

To implement the migration caused by the disk's torque F, 
we note that F = ^, and we define an extra acceleration a,„ 

2008, eq. (14)) where 



defined as a,„ = - 7- ( Cresswell & Nelson 



V is the planet's velocity, and tm = J/^ the migration time ( J is 
the angular momentum). 

In all simulations planets started on low eccentricity and low 
inclination orbits. We ran each simulation for three million years 
with a time step between 0.4 and 3 days. 



3. Two-planet case 

Figure [2] shows the evolution of two 1 planets initially on op- 
posite sides of the CZ at 3 AU. As they approach each other, the 
two planets cross a series of MMRs and are eventually trapped 
in the 7:6 resonance. The eccentricities of the two planets reach 
an equilibrium between resonant excitation and damping by the 
disk. This equilibrium eccentricity is 0.5 times the horseshoe 
semi-width x^ and damps the corotation torque to roughly 80% 
of its zero-eccentricity value. 

The planets reach a stable orbital configuration at 1.77 and 
1 .96 AU, with both planets interior to the pre-imposed CZ. Given 
their eccentricities, the innermost planet CZ is shifted to 1.95 
AU while it is 1.74 for the outer (1.96 AU) planet. In this con- 
text, the shift comes from the balance between the unchanged 
negative Lindblad torque and the attenuated positive corotation 
torque. The two-planet system stabilizes at a zero net torque 
point — even though each planet feels a nonzero torque and nei- 
ther reaches its single-planet CZ — that is shifted inward from 
the nominal CZ by > \ AU. 




0.2 



0.4 0.6 
Time [million years] 



0.8 



1.0 



Fig. 2. Simulation of the convergent migration of two 1 Mjg plan- 
ets toward the CZ at 3 AU including the Fc - e feedback illus- 
trated in the right panel of Fig. [T] 

It is clear that the planets' eccentricities — excited by planet- 
planet interactions — are the key factor in determining the 
strength of the corotation torque and the location of the eff'ective 
convergence zone. For two equal-mass planets the same qualita- 
tive behavior occurs regardless of the planet mass, although the 
evolution for higher mass planets is faster In addition, the eccen- 
tricity excitation varies both with the planet mass and which res- 
onance the planets occupy: higher eccentricities imply stronger 
damping of Fc and a closer-in effective CZ. 

3. 1 . Effect of mass ratio 

We now study the case of two unequal-mass planets. Figure |3] 
shows the outcome of a simple experiment in which a lOM® 
planet was placed at the CZ at 3 AU and a second planet was 
placed at 4 AU. The mass of the second planet was varied in a 
set of simulations from 0.1 to 10 M®. 

In Fig.[3]the outer planet remained in 3:2 resonance. In this 
situation, the final location of the planets is determined by their 
relative masses or, for this experiment, the outer planet's mass. 
The more massive the outer planet, the closer in the two planets' 
effective CZ. The inner planet's eccentricity increases for a more 
massive outer planet, causing a corresponding decrease in Fc 
and a more drastic inward shift. Given that each planet has a 
different mass and eccentricity, they have different effective CZ 
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(one per e/xj value). The magnitude of the overall inward shift, 
however, is mainly determined by the most massive planet. 

Figure [3] only shows a subset of the simulations run in this 
numerical experiment. For higher outer planet masses, the two 
planets were trapped in different MMRs, resulting in discontinu- 
ities in the equilibrium positions and mean eccentricity values. 
However, the behavior of the two planet system remained quaU- 
tatively similar. 
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Fig. 3. Outcome of a set of simulations that started with a lOMg 
planet at 3 AU and a second planet of 0.1 - 3 at 4 AU. The 
panels show the equilibrium position of the planets (top) and the 
eccentricity normalized to the horseshoe semi-width e/xs (bot- 
tom) as a function of the outer planet's mass. 



3.2. Resonance effect 

The order of MMR is also important (for an MMR of (p + q)/p, 
p is the order). Two planets in 3:2 resonance will have higher 
eccentricities than the same planets in 11:10 resonance . The 
simple explanation is that lower order resonances have more fre- 
quent conjunctions, hence stronger resonant perturbations (see 
[Murray & Dermott (2000 ) for more details). The MMR in which 
a two-planet system is captured depends on the relative migra- 
tion speed and the rate of eccentricity damping (e.g., Mustill & 
[Wyatt (2011)), both determined by the disk and torque profiles, 
and the initial positions of the planets. 

We performed a set of 100 simulations (1 Myr each) with two 
3 Mjg planets that were each initially placed randomly between 1 
and 10 AU, with the same CZ at 3 AU as before. Figure|4]shows 
that in all cases the planets are trapped in MMRs between 11:10 
and 3:2 (resonant order from 2 to 10). As expected, the excited 
eccentricity decreases for higher order resonances and leads to a 
smaller inward shift of the planetary system. The magnitude of 
the inward shift ranges from 0.2 to 1.5 AU. In two simulations 
the planets started so close to each other that they were trapped 
in co-orbital (1:1) resonance. In these cases the eccentricities re- 
mained very low, and both planets migrated to the nominal (clas- 
sical unshifted) CZ. 



4. Evolution with >2 planets 

We now turn our attention to the case of many planets in a disk. 
We ran ten simulations each with two, three, five, and ten planets 
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Fig. 4. Final semimajor axes (top) and eccentricities (bottom) of 
two 3 Me planets trapped in different MMRs. For a 3 M^ planet, 
the semi width of the horseshoe region is about 0.014. 



of 3 Mjg each. The planets were placed randomly between 1 and 
10 AU and were given low initial eccentricities and inclinations. 
As before, each simulation was evolved for 3 Myr in a static 
disk. 

For the three-planet case we find three different outcomes. 
First, the most likely outcome is that the three planets can be 
trapped in a resonant chain and migrate inward together to a zero 
net torque zone. This zone is typically between 2 and 2.5 AU. 
The eccentricities of the three planets are not the same: the mid- 
dle planet is generally more excited. In the second most likely 
outcome, two planets are trapped in MMR and shifted inward, 
but a third, outer planet is too far away to be trapped in a reso- 
nant chain. The outer planet remains at the nominal CZ at 3 AU, 
while the inner planets end up at a new shifted zero torque zone. 
Finally, in some cases collisions occur between two planets and 
the system reverts to the two planet, unequal-mass case seen in 
Fig. I 

For the case of five or ten planets the situation is more 
complex. The five-planet systems become trapped in chains 
of MMRs and migrate inward to an equilibrium position with 
zero net torque. However, perturbations between planets add a 
stochastic element to the planets' behavior. Even the most sta- 
ble systems include periods during which the planets all drift 
radially in the same direction. These periods are triggered by 
the breaking of a resonance between two planets that propagates 
as a perturbation through the whole system. The amplitude and 
frequency of these chaotic periods vary from case to case. For 
example, in the top simulation from Fig. |5] the resonant chain 
undergoes several small hiccups but the drift is small compared 
with the typical inter-planetary spacing. In contrast, the pertur- 
bations in the simulation in the bottom panel of Fig.|5]are much 
stronger For instance, we consider the interval between roughly 
1.1 and 1.3 Myr in the bottom case. At 1.12 Myr the two outer 
planets become trapped in 4:3 MMR. They migrate inward due 
to the resonant eccentricity increase. This perturbation propa- 
gates to the inner system, and the eccentricity increase causes the 
planets to migrate inward all together. Five thousand years later, 
the outer planets leave and then re-enter the 4:3 MMR, again 
perturbing the system. Finally, at 1.13 Myr the outer two planets 
leave the 4:3 MMR for good. As they leave the MMR, the two 
planets act like isolated bodies and migrate toward the nominal 
convergence zone with near zero eccentricities. Without the net 
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negative torque supplied by the outer two planets, the inner plan- 
ets reactively migrate outward toward an effective three-planet- 
system zero torque zone. This is the cause of the abrupt outward 
migration of all five planets. However, the two outer planets 
quickly find themselves in the 5:4 MMR. Over the subsequent 
0.15 Myr interval, they are periodically trapped in 5:4 MMR but 
outward migration continues because most of the time they are 
not in MMR and their eccentricity is low. The system-wide out- 
ward migration stops at 1.335 Myr when the outer planets cross 
over the 5:4 MMR and are trapped in 6:5 MMR. This configura- 
tion stabilize the system, excites the outer planets' eccentricities, 
and causes the whole system to migrate inward again, which sig- 
nals the end of this particular meta perturbation. 

The rest of the evolution is composed of the same kinds of 
perturbations. Perturbations come from planets entering or leav- 
ing resonances, then propagate between planets. When leaving 
a resonance, planets' eccentricities usually drop, causing out- 
ward migration. Conversely, being trapped in resonance excites 
a sustained eccentricity that leads to inward migration. Via these 
perturbations, entire systems of planets undergo modest-scale 
chaotic migration due to the difficulty of maintaining resonant 
chains over long times. Each of the five planet systems we sim- 
ulated remained stable, in that no collisions occurred, but the 
amplitude of the chaotic migration varied from case to case; the 
two examples from Fig. [5] show the extremes. The ten-planet 
simulations were even more chaotic, and collisions did occur 

The most important factor in determining the ampUtude of 
chaotic oscillations is the order of the resonances. Low-order 
resonances maintain higher eccentricities and are less stable be- 
cause they are sensitive to eccentricity variations. On the other 
hand, high-order resonances maintain lower eccentricities, and 
eccentricity perturbations are less important. For instance, in 
the system in the top panel of Fig. |5] perturbations are quickly 
smoothed out, whereas in the system in the bottom panel the fre- 
quency of perturbations is high enough that the system does not 
tend to a stably evolving configuration. Therefore, in this con- 
text a more compact system of planets is more stable than an 
extended system, which is the opposite of the pure gravitational 
situation ( |Marchal & Bozis|1982j . 




1.0 1.5 2.0 

Time [million years] 



3.0 





0.0 



0.5 



xxi rs 2J5 

Time [million years 



yAriM^V^^ References 



5. Discussion 

We have shown that planets do not actually converge in con- 
vergence zones (CZs). Instead, embryos rapidly migrate toward 
the CZ and are trapped in chains of MMRs. This causes their 
eccentricities to increase and remain high enough to attenuate 
the corotation torque. The zero-torque equilibrium point of the 
resonant chain of planets is determined by the sum of each in- 
dividual planet torque (sum of attenuated corotation torque and 
unattenuated differential Linbdlad torque). In practice, the effec- 
tive zero net torque zone is shifted inward and is most strongly 
determined by the CZ of the most massive planet in the resonant 
chain. It is not a true CZ because each planet feels a different CZ 
(depending on its eccentricity). 

The inward shift exists because planets' eccentricities are 
sustained by resonant perturbations. The amplitude of a planet's 
eccentricity is a result of the competition between resonant exci- 
tation and eccentricity damping by the disk. For sufficiently high 
sustained eccentricities, an entire resonant chain of planets can 
migrate all the way in to the inner edge of the disk. Changing the 
disk properties may thus change the typical sustained eccentric- 
ities by affecting the eccentricity damping timescale. However, 
changing the disk properties affects other properties of the sys- 
tem such as the torque profile. In changing the properties of the 
disk, it is not straightforward to assess the effect on the plan- 
ets' evolution, as planets may be trapped in different resonances 
with different eccentricity excitations and even different stability 
criteria. 

The CZ depends on the disk parameters such as the viscos- 



ity, temperature and surface density profiles (e.g. Paardekooper 
|et al.|20lT] l. Here we have used a disk profile that is motivated 
by complex models but is nonetheless artificial. Even if the re- 
sults shown use a particular model, they are robust when we vary 
the torque shape in function of the orbital distance, as long as a 
CZ zone exists. In a more realistic disk we would expect a few 
differences. First, there may be multiple CZs in the same disk 
owing to different effects ( Lyra et al.|2010[ Hasegawa & Pudritz 



201 1 ). Second, mass-dependent convergence zones may exist in 
the outer part of the disk, where this mechanism should be less 
efficient because planets do not migrate to the same location in 
the disk and may not produce the resonant chains essential for 
our mechanism to occur (Cossou et al, in preparation). Third, as 
the disk dissipates, so does the torque profile and the location 
of CZs ( [Lyra et al.|[20Tm [Horn et al.|20T2] ). Finally, turbulence 
is thought to be common in real protoplanetary disks ( |Armitage| 
201 1 1. Although turbulence does not affect the long-term evolu- 
tion of a single planet in a radiative disk ( jPierens et al. 2012 1, we 
expect it to modify resonant capture and eccentricity evolution 
(see |Pierens et al.|201 \) . 
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Fig. 5. Two example simulations with 5 planets, one that is quite 
stable even if short resonance breaking occurs (top panel) and 
one that exhibits sustained chaotic behavior (bottom). 
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